1 Introduction

This file was last update on 2020-10-02

2 Sampling effort

No. of visits to each transect/site.

abn%>%
  group_by(Study, Taxa, Transect)%>%
  summarise(n = length(unique(Date)))%>%
  spread(Transect, n)
Study Taxa T1A T1B T2 T3 T4 T5 T6
1998 - 2000 Birds 19 18 7 10 4 1 1
1998 - 2000 Butterflies 14 7 4 9 4 1 1
2016 - 2017 Birds 23 23 22 19 16 22 17
2016 - 2017 Butterflies 23 20 20 22 18 21 15

No. of individuals sampled

abn%>%
  group_by(Study, Taxa, Transect)%>%
  summarise(n = sum(Abundance))%>%
  spread(Transect, n)
Study Taxa T1A T1B T2 T3 T4 T5 T6
1998 - 2000 Birds 255 165 103 213 236 7 28
1998 - 2000 Butterflies 93 59 12 180 157 10 4
2016 - 2017 Birds 480 393 358 289 165 203 133
2016 - 2017 Butterflies 432 188 207 652 187 224 124

Uneveness in sampling effort across sites in studies may bias estimates of change.

2.1 Rarefying across sites

library(vegan)

# input for vegan

veg_birds <- abn%>%
  filter(Taxa == "Birds")%>%
  group_by(Study, Transect, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = paste(Study, "_", Transect))%>%
  group_by(sample)%>%
  #mutate(N = sum(n, na.rm = T))%>%
  #filter(N > 0)%>%#removing samples where no species were detected
  dplyr::select(sample, Specific.Epithet, n)%>%
  ungroup()%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

veg_butterfly <- abn%>%
  filter(Taxa == "Butterflies")%>%
  group_by(Study, Transect, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = paste(Study, "_", Transect))%>%
  group_by(sample)%>%
  #mutate(N = sum(n, na.rm = T))%>%
  #filter(N > 0)%>%#removing samples where no species were detected
  dplyr::select(sample, Specific.Epithet, n)%>%
  ungroup()%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")
## Species accumulation curves

specaccum(veg_butterfly)
## Species Accumulation Curve
## Accumulation method: exact
## Call: specaccum(comm = veg_butterfly) 
## 
##                                                                        
## Sites     1.00000  2.00000  3.00000  4.00000  5.00000  6.00000  7.00000
## Richness 26.85714 39.51648 47.34615 52.96004 57.33017 60.91275 63.96037
## sd       12.83331 11.06085  9.33188  8.00650  7.02146  6.19857  5.49616
##                                                                  
## Sites     8.00000  9.00000 10.00000 11.00000 12.00000 13.00000 14
## Richness 66.62438 68.99800 71.13986 73.08791 74.86813 76.50000 78
## sd        4.88039  4.31448  3.76625  3.21657  2.59888  1.76271  0
rareslope(veg_butterfly, sample = 200)
## 1998 - 2000 _ T1A 1998 - 2000 _ T1B  1998 - 2000 _ T2  1998 - 2000 _ T3 
##        0.00000000        0.00000000        0.00000000        0.00000000 
##  1998 - 2000 _ T4  1998 - 2000 _ T5  1998 - 2000 _ T6 2016 - 2017 _ T1A 
##        0.00000000        0.00000000        0.00000000        0.04559341 
## 2016 - 2017 _ T1B  2016 - 2017 _ T2  2016 - 2017 _ T3  2016 - 2017 _ T4 
##        0.00000000        0.05079708        0.04937664        0.00000000 
##  2016 - 2017 _ T5  2016 - 2017 _ T6 
##        0.05173746        0.00000000
rarefy(veg_butterfly, sample = 200)
## 1998 - 2000 _ T1A 1998 - 2000 _ T1B  1998 - 2000 _ T2  1998 - 2000 _ T3 
##          27.00000          22.00000           8.00000          37.00000 
##  1998 - 2000 _ T4  1998 - 2000 _ T5  1998 - 2000 _ T6 2016 - 2017 _ T1A 
##          24.00000           5.00000           3.00000          34.98930 
## 2016 - 2017 _ T1B  2016 - 2017 _ T2  2016 - 2017 _ T3  2016 - 2017 _ T4 
##          32.00000          35.65386          31.83323          31.00000 
##  2016 - 2017 _ T5  2016 - 2017 _ T6 
##          37.84749          27.00000 
## attr(,"Subsample")
## [1] 200
rarebutts <- rarecurve(veg_butterfly)

3 Change in diversity across years

3.1 Species richness

3.1.1 Across the two studies

abn%>%
  group_by(Taxa, Study)%>%
  summarise(richness = length(unique(Specific.Epithet)),
            N = sum(Abundance, na.rm = T))
Taxa Study richness N
Birds 1998 - 2000 70 1007
Birds 2016 - 2017 105 2021
Butterflies 1998 - 2000 45 515
Butterflies 2016 - 2017 66 2014

More species of both taxa were encountered in current study compared to the previous study. Difference in sampling effort is apparent accross two studies.

3.1.2 Across sites in the two studies

No. of species encountered across sites sampled in the two studies.

abn%>%
  left_join(transects, by = c("Transect"))%>%
  group_by(Taxa, Study, Transect)%>%
  summarise(richness = length(unique(Specific.Epithet)))%>%
  spread(Transect, richness)
Taxa Study T1A T1B T2 T3 T4 T5 T6
Birds 1998 - 2000 47 37 28 36 29 4 11
Birds 2016 - 2017 66 59 53 47 31 38 28
Butterflies 1998 - 2000 27 22 8 37 24 5 3
Butterflies 2016 - 2017 41 32 36 44 31 39 27

Fewer species were encountered in sites with lower effort in the previous study, i.e., T5 and T6.

3.2 Shannon diversity

# Calculating diversity in each site over years

div_year <- abn%>%
  group_by(Study, Taxa, Transect,  Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  spread(Specific.Epithet, n, fill = c(0))%>%
  group_by(Study, Taxa, Transect)%>%
  nest()%>%
  mutate(H = map_dbl(data, ~vegetarian::d(.[-1], lev = "alpha", q = 1)))%>%
  dplyr::select(-data)

# Summary

div_year%>%
  group_by(Taxa, Study)%>%
  skimr::skim(H)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 15.39 7.79
Birds 2016 - 2017 24.88 5.74
Butterflies 1998 - 2000 11.06 6.57
Butterflies 2016 - 2017 20.00 3.01 
# t test

div_year%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(H ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa 1998-2001 2016-2017 parameter statistic p.value
Birds 15.39287 24.87583 11.031229 -2.593073 0.0249503
Butterflies 11.05543 19.99721 8.410531 -3.273453 0.0105504 
# plot

div_year%>%
  ggplot(aes(Taxa, H, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Effective number of species")+
  scale_fill_brewer(palette = "Greys")

ggsave(last_plot(), filename = "./Figures and Tables/figure2.png", height = 6, width = 8, dpi = 300)

The first order diversity of both taxa has increased significantly across the two surveys.

4 Change in composition of species accross years

# sample info

samp_birds <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  filter(Taxa == "Birds")%>%
  dplyr::select(Study, Taxa, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))

samp_butterfly <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  filter(Taxa == "Butterflies")%>%
  dplyr::select(Study, Taxa, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))
# Similarity matrix by studies

sim_birds <- abn%>%
  filter(Taxa == "Birds")%>%
  group_by(Study, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = Study)%>%
  ungroup()%>%
  dplyr::select(sample, Specific.Epithet, n)%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

sim_butterfly <- abn%>%
  filter(Taxa == "Butterflies")%>%
  group_by(Study, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = Study)%>%
  ungroup()%>%
  dplyr::select(sample, Specific.Epithet, n)%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

# bray - curtis

data_frame(
  Taxa = c("Birds", "Butterflies"),
  Dissimilarity = c(as.numeric(vegdist(sim_birds, method = "bray")),
                    as.numeric(vegdist(sim_butterfly, method = "bray")))
)
Taxa Dissimilarity
Birds 0.5495376
Butterflies 0.6710162 
# ordination

NMDS_birds <- metaMDS(veg_birds, distance = "bray", k = 2, model = "local")

NMDS_butterfly <- metaMDS(veg_butterfly, distance = "bray", trymax = 1000, k = 2, model = "local")


# Extracting scores and adding habitat vars

ord_birds <- data.frame(scores(NMDS_birds))%>%
  rownames_to_column("sample")%>%
  left_join(samp_birds, "sample")

ord_butterfly <- data.frame(scores(NMDS_butterfly))%>%
  rownames_to_column("sample")%>%
  left_join(samp_butterfly, "sample")

ord <- bind_rows(ord_birds, ord_butterfly)
#Testing change in composition of birds

adonis(veg_birds ~ Study, data = samp_birds)
## 
## Call:
## adonis(formula = veg_birds ~ Study, data = samp_birds) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## Study      1   0.77547 0.77547  3.9632 0.24827  0.002 **
## Residuals 12   2.34798 0.19567         0.75173          
## Total     13   3.12345                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Testing change in composition of butterflies

adonis(veg_butterfly ~ Study, data = samp_butterfly)
## 
## Call:
## adonis(formula = veg_butterfly ~ Study, data = samp_butterfly) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Study      1    0.9191 0.91908   4.054 0.25253  0.001 ***
## Residuals 12    2.7205 0.22671         0.74747           
## Total     13    3.6396                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# plotting

ord%>%
  ggplot(aes(NMDS1, NMDS2, col = Study, fill = Study))+
  stat_ellipse(level = 0.95, geom = "polygon", linetype = "dashed", size = 1, alpha = 0.3)+
  geom_point(aes(shape = Habitat), size = 3)+
  guides(shape = guide_legend(nrow = 2, byrow = T))+
  scale_color_brewer(palette = "Set1")+
  scale_fill_brewer(palette = "Set1")+
  facet_wrap(~Taxa)+
  theme(legend.box = "vertical")

ggsave(last_plot(), filename = "./Figures and Tables/figure3.png", height = 6, width = 8, dpi = 300)

Both bird and butterfly communities compositions are significantly different compared to the previous study.

5 Change in niche width of community

5.1 Trophic Guilds

# compiling trophic guild data

host <- read.csv("./Data/Butterflies_host plant_280820.csv") #butterfly host plant data

diet <- read.csv("./Data/Birds_diet.csv") #bird diet data

trophic <- host%>%
  full_join(diet)%>%
  dplyr::select(Specific.Epithet, Guild)%>%
  distinct()

# Calculating diversity of guilds

guilds <- abn%>%
  left_join(trophic, by = c("Specific.Epithet"))%>%
  filter(!Guild == "")%>%
  group_by(Study, Taxa, Transect, Guild,  Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  spread(Specific.Epithet, n, fill = c(0))%>%
  group_by(Study, Taxa, Transect, Guild)%>%
  nest()%>%
  mutate(H = map_dbl(data, ~vegetarian::d(.[-1], lev = "alpha", q = 1)))%>%
  dplyr::select(-data)

# summary

guilds%>%
  group_by(Taxa, Guild, Study)%>%
  skimr::skim(H)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Guild, Study, mean, sd)

Variable type: numeric

Taxa Guild Study mean sd
Birds Carnivore 1998 - 2000 3.01 1.14
Birds Carnivore 2016 - 2017 4.67 2.62
Birds Frugivore 1998 - 2000 3.37 1.51
Birds Frugivore 2016 - 2017 4.21 0.60
Birds Grainivore 1998 - 2000 2.44 1.43
Birds Grainivore 2016 - 2017 2.83 1.09
Birds Insectivore 1998 - 2000 7.74 3.06
Birds Insectivore 2016 - 2017 13.58 4.82
Birds Omnivore 1998 - 2000 3.90 1.78
Birds Omnivore 2016 - 2017 5.88 1.15
Butterflies Generalist 1998 - 2000 2.79 1.28
Butterflies Generalist 2016 - 2017 5.62 1.91
Butterflies Grass Specialist 1998 - 2000 1.87 0.57
Butterflies Grass Specialist 2016 - 2017 2.98 0.74
Butterflies Herb Specialist 1998 - 2000 2.73 0.93
Butterflies Herb Specialist 2016 - 2017 4.09 0.86
Butterflies Liana Specialist 1998 - 2000 1.00 0.00
Butterflies Shrub Specialist 1998 - 2000 3.54 1.71
Butterflies Shrub Specialist 2016 - 2017 4.78 1.50
Butterflies Tree Specialist 1998 - 2000 4.63 3.02
Butterflies Tree Specialist 2016 - 2017 6.88 1.63 
# t - test

guilds%>%
  filter(Guild != "Liana Specialist")%>% #removing due to lack of obs
  complete(nesting(Study, Taxa), Transect, Guild, fill = list(H = 0))%>%
  group_by(Taxa, Guild)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(H ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa Guild 1998-2001 2016-2017 parameter statistic p.value
Birds Carnivore 3.013676 3.334886 40.52328 -0.6191168 0.5393047
Birds Frugivore 3.367593 4.211182 62.64226 -3.8933034 0.0002426
Birds Grainivore 2.440818 2.833390 89.94433 -1.6340146 0.1057514
Birds Insectivore 6.631103 13.578874 88.97385 -7.9838426 0.0000000
Birds Omnivore 3.901475 5.876942 82.45283 -6.9810904 0.0000000
Butterflies Generalist 1.990058 5.622563 77.56402 -9.7435153 0.0000000
Butterflies Grass Specialist 1.335180 2.975740 57.94419 -8.6211871 0.0000000
Butterflies Herb Specialist 1.560185 4.092249 35.98571 -8.2392106 0.0000000
Butterflies Shrub Specialist 3.542509 4.780687 94.34174 -4.0786608 0.0000946
Butterflies Tree Specialist 4.631978 6.879746 73.68039 -4.8985262 0.0000056 
# plot
guilds%>%
  ggplot(aes(reorder(Guild, H), H, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5, position = position_dodge(preserve = "single"))+
  labs(y = "Effective number of species", x = "Trophic Guild")+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  scale_fill_brewer(palette = "Greys")+
  facet_wrap(~Taxa, scales = "free_x", ncol = 1)

ggsave(last_plot(), filename = "./Figures and Tables/figure4.png", height = 9, width = 8, dpi = 300)

Carnivorous, Grainivorous birds have reduced in the comunity, compared to insectivorous birds which have increased. Similarly, Shrub specialist and tree specialist butterflies have reduced where grass specialists and herb specialist have increased. Suprsingly, generalist species of both taxa appear to be unchanged.

6 Change in niche specialisation

-Can we incorporate abunndance into community specialisation index?

6.1 Trophic niche

# trophic specialisation based on Juliard et al. 2006

butterflies_TSI <- host%>%
  mutate(H = length(unique(Hostplant.Family)))%>%
  group_by(Specific.Epithet)%>%
  mutate(h = length(unique(Hostplant.Family)))%>%
  ungroup()%>%
  mutate(TSI = sqrt((H/h)-1))

birds_TSI <- diet%>%
  mutate(H = length(unique(Prey)))%>%
  group_by(Specific.Epithet)%>%
  mutate(h = length(unique(Prey)))%>%
  ungroup()%>%
  mutate(TSI = sqrt((H/h)-1))

CSI.T <- abn%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
  group_by(Study, Taxa, Transect)%>%
  distinct()%>%
  left_join(butterflies_TSI, by = c("Specific.Epithet"))%>%
  left_join(birds_TSI, by = c("Specific.Epithet"))%>%
  mutate(TSI = ifelse(is.na(TSI.x), TSI.y, TSI.x))%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet, TSI)%>%
  summarise(CSI.T = mean(TSI, na.rm = T))
# summary
CSI.T%>%
  group_by(Taxa, Study)%>%
  skimr::skim(CSI.T)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 2.59 0.16
Birds 2016 - 2017 2.53 0.09
Butterflies 1998 - 2000 5.34 0.39
Butterflies 2016 - 2017 5.07 0.24 
# t - test

CSI.T%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(CSI.T ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa 1998-2001 2016-2017 parameter statistic p.value
Birds 2.588009 2.531490 9.462567 0.8103492 0.4376441
Butterflies 5.338563 5.065191 10.056310 1.5727722 0.1466765 
CSI.T%>%
  ggplot(aes(Taxa, CSI.T, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Community specialisation index \n (Trophic)")+
  scale_fill_brewer(palette = "Greys")+
  scale_y_sqrt()

Niche width of the community did not change over the two studies.

6.2 Habitat Niche

HSI <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  group_by(Study, Taxa, Habitat, Specific.Epithet)%>%
  # calculating abundance of each species in each habitat type
  summarise(n = sum(Abundance, na.rm = T))%>%
  group_by(Study, Taxa)%>%
  mutate(N = sum(n), # total number of indviduals encountered 
         rel.prop = n/N)%>% # relative proportion of species
  group_by(Study, Specific.Epithet)%>%
  summarise(HSI = sd(rel.prop)/mean(rel.prop))%>%
  replace_na(replace = list(HSI = 0))

CSI.H <- abn%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
  group_by(Study, Taxa, Transect)%>%
  distinct()%>%
  left_join(HSI, by = c("Study", "Specific.Epithet"))%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet, HSI)%>%
  summarise(CSI.H = mean(HSI, na.rm = T))
# summary
CSI.H%>%
  group_by(Taxa, Study)%>%
  skimr::skim(CSI.H)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 0.54 0.10
Birds 2016 - 2017 0.49 0.02
Butterflies 1998 - 2000 0.61 0.09
Butterflies 2016 - 2017 0.57 0.03 
# t - test

CSI.H%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(CSI.H ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()
Taxa 1998-2001 2016-2017 parameter statistic p.value
Birds 0.5363977 0.4942803 6.719090 1.068760 0.3220632
Butterflies 0.6110527 0.5692785 7.782004 1.210166 0.2616879 
CSI.H%>%
  ggplot(aes(Taxa, CSI.H, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Community specialisation index \n (Habitat)")+
  scale_fill_brewer(palette = "Greys")

Habitat niche of the community did not change over the study periods.

7 Seasonal changes in relative abundance

season <- abn%>%
  mutate(Month = month(dmy(Date), label = T))%>%
  group_by(Study, Month, Taxa)%>%
  summarise(n = sum(Abundance, na.rm=T))%>%
  group_by(Study, Taxa)%>%
  mutate(N = sum(n),
         rel.prop = n/N)

season%>%
  ggplot(aes(Month, rel.prop, col = Taxa, shape = Taxa))+
  geom_point(size = 3)+
  geom_path(aes(group = Taxa), size = 1, linetype = "dashed")+
  scale_color_brewer(palette = "Set1")+
  facet_wrap(~Study, ncol = 1)

8 Site Map

library(raster)
library(sf)

# Raw point data

points <- read.csv("./Data/sites.csv")%>%
  drop_na()%>%
  dplyr::select(-Location)

# Converting to SF points

points_sf <- st_as_sf(points, coords = c('Lon', 'Lat'))

st_crs(points_sf) <- 4326

# Making transect lines

lines <- points_sf%>%
  mutate(Transect = as.factor(Transect))%>%
  group_by(Transect)%>%
  summarise()%>%
  st_cast("LINESTRING")
# getting study extent

study_ext <- st_bbox(points_sf)

ext <- st_as_sfc(study_ext)

# getting osm map

library(osmdata)

# Tamhini WLS

tamhini_raw <- opq(bbox = study_ext, timeout = 100)%>%
  add_osm_feature(key = 'boundary', value = 'protected_area')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
tamhini <- tamhini_raw$osm_multipolygons

st_crs(tamhini) <- 4326

# Roads

tamhini_area <- st_bbox(tamhini)

roads_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'highway')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
roads <- roads_raw$osm_lines

st_crs(roads) <- 4326

# Landuse

landuse_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'natural')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
forest <- landuse_raw$osm_polygons%>%
  filter(natural == "wood")

water <- landuse_raw$osm_multipolygons

water_small <- landuse_raw$osm_polygons%>%
  filter(natural == "water")

st_crs(forest) <- 4326

# waterways

waterway_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'waterway')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
waterways <- waterway_raw$osm_lines

st_crs(waterways) <- 4326

# elevation 

library(elevatr)

tam_elev <- elevatr::get_elev_raster(locations = tamhini, z = 8)

tam_contours <- raster::rasterToContour(tam_elev, levels = seq(0, 1000, 100))
# making final map

library(gridExtra)

fig1 <- grid.arrange(fig1a, fig1b, fig1c, legend, layout_matrix = rbind(c(1, 1, 1),
                                                          c(1, 1, 1),
                                                          c(2, 3, 4)))

ggsave(fig1, filename = "./Figures and Tables/figure1.png", height = 8, width = 8)